subroutine write_hdf5_data(isnap,iframe,itype)
  use logs
  use my_mpi
  use parameter_module
  use camera_module
  use grid
  use smoothing_variables
  use myfiles
  use hdf5_wrapper
  use header
  use units_and_constants
  implicit none
  integer*4 :: i,isnap,iframe,itype,ivar
  integer*4 :: ierr,iPE,ifile,file_handle,ntotal,ioffset,nremain
  integer*4, parameter :: rank = 7
  integer*4, dimension(rank) :: count  
  integer*4, dimension(rank) :: offset
  integer*4, dimension(rank) :: init_size

  real*4, allocatable :: temp2D(:,:),out2D(:,:)
  real*4, allocatable :: temp3D(:,:,:)

  character(len=256) :: framedir,speciesdir
  character(len=256) :: outpath,outfile,title,dir
  character(len=16) :: ctype

  write(logunit,*)'  -writing HDF5 data...'
  call flush(logunit)

  call mpi_barrier(mpi_comm_world,ierr)

  ! Set the directory name
  write(framedir,'("Grid_",I3.3)')iframe
  write(speciesdir,'("PartType",I1)')itype
  outpath = trim(params%outputdir)//'/'//trim(framedir)//'/'//trim(speciesdir)   

  ! Generate this directory 
  call system('mkdir -p '//trim(outpath))

  ! Set the output filename
  outfile  = trim(outpath)//'/Grids.hdf5'
  write(logunit,*)'   -outfile: '//trim(adjustl(outfile))

  ! Add optional title to vtk file
  title = 'Made with GridHDF5 :: Rob Crain, ICC, Durham'

  ! Create HDF5 file on root process, then close it so that all procs 
  ! then open it and have access to the file handle
  if(MyPE==0)then
     call hdf5_create_file(file_handle,trim(outfile))

     ! Write metadata
     call hdf5_create_group(file_handle,'MetaData')
#ifdef GADGET2
     call hdf5_write_attribute(file_handle,'MetaData/Flag_Gadget2',1)
#else
     call hdf5_write_attribute(file_handle,'MetaData/Flag_Gadget2',0)
     call hdf5_write_attribute(file_handle,'MetaData/Time_GYR',Time_GYR)
#endif
     call hdf5_write_attribute(file_handle,'MetaData/Title',trim(title))
     call hdf5_write_attribute(file_handle,'MetaData/isnap',isnap)
     call hdf5_write_attribute(file_handle,'MetaData/iframe',iframe)
     call hdf5_write_attribute(file_handle,'MetaData/Redshift',redshift)
     call hdf5_write_attribute(file_handle,'MetaData/Centre_in_CodeUnits',(/x0,y0,z0/))
#ifdef VELOCITY
     call hdf5_write_attribute(file_handle,'MetaData/Velocity_in_CodeUnits',(/vx0,vy0,vz0/))
#endif
     call hdf5_write_attribute(file_handle,'MetaData/L',L(:))
#ifdef USE_EIGENVECTORS
     call hdf5_write_attribute(file_handle,'MetaData/Flag_EigenVectors',1)
     call hdf5_write_attribute(file_handle,'MetaData/X-Prime',xprime)
     call hdf5_write_attribute(file_handle,'MetaData/Y-Prime',yprime)
     call hdf5_write_attribute(file_handle,'MetaData/Z-Prime',zprime)
#else
     call hdf5_write_attribute(file_handle,'MetaData/Flag_EigenVectors',0)
     call hdf5_write_attribute(file_handle,'MetaData/Theta',theta)
     call hdf5_write_attribute(file_handle,'MetaData/Phi',phi)
     call hdf5_write_attribute(file_handle,'MetaData/Psi',psi)
#endif
#ifdef CUBE
     call hdf5_write_attribute(file_handle,'MetaData/Ngrid',(/params%nx,params%ny,params%nz/))
#else
     call hdf5_write_attribute(file_handle,'MetaData/Ngrid',(/params%nx,params%nz/))
#endif
     call hdf5_write_attribute(file_handle,'MetaData/ComovingDistanceScale_Mpc',Dcom)
     call hdf5_write_attribute(file_handle,'MetaData/ComovingAngularScale_kpc_per_arcsec',comoving_angular_scale)
     call hdf5_write_attribute(file_handle,'MetaData/PhysicalAngularScale_kpc_per_arcsec',physical_angular_scale)

     ! Write Header data
     call write_header(file_handle)

     ! Close file
     call hdf5_close_file(file_handle)
  endif

  ! Barrier here so non-root processes don't try to open the output file 
  ! before the root creates it
  call mpi_barrier(mpi_comm_world,ierr)

  ! Cycle over the cores and write the dataset independently on each CPU.
  do iPE=0,NumPEs-1

     if(iPE==MyPE) then

        ! Select hyperslab in the file.
        NZLocal  = params%nz / NumPEs
        NRemain  = NZLocal*(NumPEs-1)
        NZMaster = params%nz - NRemain
        if(MyPE==0) NZLocal = NZMaster

#ifdef CUBE
        offset(:) = 1
        count(:)  = 1
        offset(3) = NZMaster + (MyPE-1)*NZLocal + 1
        count(1)  = params%nx
        count(2)  = params%ny
        count(3)  = nzlocal
        init_size = (/params%nx,params%ny,params%nz,1,1,1,1/)
#else
        offset(:) = 1
        count(:)  = 1
        offset(2) = NZMaster + (MyPE-1)*NZLocal + 1
        count(1)  = params%nx
        count(2)  = NZLocal
        init_size = (/params%nx,params%nz,1,1,1,1,1/)
#endif

        call hdf5_open_file(file_handle,trim(adjustl(OutFile)),readonly=.false.)

        ! Write weighting grid
        dir = 'Grids/'//trim(adjustl(speciesdir))//'/Weight'
        call hdf5_write_data(file_handle,trim(adjustl(dir)),Value,&
             start=offset,count=count,initial_size=init_size)
        if(MyPE==0)then
           call hdf5_write_attribute(file_handle,trim(dir)//'/Minimum',ValueMin)
           call hdf5_write_attribute(file_handle,trim(dir)//'/Maximum',ValueMax)
           call hdf5_write_attribute(file_handle,trim(dir)//'/Total'  ,ValueSum)
        endif

        do ivar=1,nvar_by_type(itype)

           if(trim(var(itype,ivar))=='PartType0/Weight') cycle	
           if(trim(var(itype,ivar))=='PartType1/Weight') cycle	
           if(trim(var(itype,ivar))=='PartType2/Weight') cycle	
           if(trim(var(itype,ivar))=='PartType3/Weight') cycle	
           if(trim(var(itype,ivar))=='PartType5/Weight') cycle	
           if(trim(var(itype,ivar))=='PartType6/Weight') cycle	
           if(trim(var(itype,ivar))=='PartType7/Weight') cycle

           dir = 'Grids/'//trim(adjustl(var(itype,ivar)))
#ifdef CUBE
           allocate(temp3D(params%nx,params%ny,NZLocal))
           temp3D = ValueQuantity(ivar,:,:,:)
           call hdf5_write_data(file_handle,trim(adjustl(dir)),temp3D,&
                start=offset,count=count,initial_size=init_size)
           deallocate(temp3D)
#else
           allocate(temp2D(params%nx,NZLocal))
           temp2D = ValueQuantity(ivar,:,:)
           call hdf5_write_data(file_handle,trim(adjustl(dir)),temp2D,&
                start=offset,count=count,initial_size=init_size)
           deallocate(temp2D)

#endif
           if(MyPE==0)then
              call hdf5_write_attribute(file_handle,trim(dir)//'/Minimum',ValueQuantityMin(ivar))
              call hdf5_write_attribute(file_handle,trim(dir)//'/Maximum',ValueQuantityMax(ivar))
              call hdf5_write_attribute(file_handle,trim(dir)//'/Total'  ,ValueQuantitySum(ivar))
           endif

        enddo

        call hdf5_close_file(file_handle)        

     endif

     call mpi_barrier(mpi_comm_world,ierr)
  enddo

  write(logunit,*)'   -finished writing'
  call flush(logunit)

end subroutine write_hdf5_data
